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Elastic-plastic  behavior  of  cyclotrimethylene  trinitramine  single  crystals 
under  spherical  indentation:  Modeling  and  simulation 
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(Received  5  January  2012;  accepted  13  February  2012;  published  online  21  March  2012) 

A  nonlinear  anisotropic  elastic-plastic  model  is  developed  for  single  crystals  of  the  energetic 
material  cyclotrimethylene  trinitramine  (RDX).  Numerical  simulations  of  spherical  indentation 
on  oriented  single  crystals  are  performed;  predictions  are  compared  with  experimental  data  and 
observations  from  the  literature.  Model  predictions  for  elastic  response  and  initial  yield  using 
elastic  constants  obtained  from  resonant  ultrasound  spectroscopy  agree  with  experimental  data; 
predicted  forces  using  constants  obtained  from  Brillouin  scattering  tend  to  exceed  experimental 
data.  Influences  of  elastic  anisotropy  and  elastic  nonlinearity  are  significant.  Predicted  slip 
system  activity  is  in  reasonable  agreement  with  that  deduced  from  experimental  surface  profiles 
when  a  uniform  strength  of  G/ 20  is  assigned  to  all  six  slip  systems,  with  G  an  effective  elastic 
shear  modulus.  Predicted  indentation  forces  in  the  post-yield  regime  exceed  those  observed  in 
experiments,  suggesting  that  surface  and  possibly  subsurface  fractures  may  contribute  to  a  loss 
of  stiffness  in  experiments  at  larger  indentation  depths.  [http://dx.doi.org/10.1063/L3695392] 


I.  INTRODUCTION 

Defects  in  energetic  materials  are  thought  to  affect  their 
initiation  sensitivity.  Stresses  concentrate  in  the  vicinity  of 
cracks,  pores,  or  lattice  defects,  which  in  turn  can  affect  ini¬ 
tiation  of  reactions  associated  with  burning  or  detonation.  In 
single  crystals,  availability  of  slip  systems  associated  with 
mobile  dislocations  may  lower  peak  stresses  and  decrease 
sensitivity  to  shock  initiation. 1 

The  focus  of  the  present  work  is  the  mechanical  behav¬ 
ior  of  the  energetic  material  cyclotrimethylene  trinitramine 
(C3H6N606),  referred  to  as  RDX  (Research  Development 
eXplosive).  Single  crystals  of  RDX  belong  to  an  orthorhom¬ 
bic  space  group  with  eight  molecules  per  unit  cell.  Disloca¬ 
tions  in  RDX  have  been  characterized  using  etch  pit2  and  x 
ray  topographic3'4  techniques.  Likely  slip  systems  in  RDX 
have  been  suggested  from  analysis  of  anisotropic  hardness 
profiles5  and  indentation  experiments.6'7  The  latter  experi¬ 
ments6'7  also  provide  an  estimate  of  the  critical  resolved 
shear  stress  associated  with  slip  initiation,  thought  to  be  on 
the  order  of  the  theoretical  strength  (i.e.,  «  G/10-G/20,  with 
G  a  representative  elastic  shear  modulus),  which  corresponds 
to  homogeneous  dislocation  nucleation.  Inelastic  behavior  of 
RDX  crystals  has  also  been  probed  using  shock  experi¬ 
ments1'8  and  molecular  dynamics  simulations.8'9 

Continuum  crystal  plasticity  theory  permits  predictive 
mesoscale  modeling  of  materials’  behavior  at  length  scales 
larger  than  that  feasible  using  molecular  models,  but  with 
greater  resolution  than  that  afforded  by  macroscopic  elastic- 
plastic  models  that  omit  anisotropy  and  slip  system  activity. 
Grain  interactions  can  be  studied  in  direct  numerical  simula¬ 
tions  via  finite  element  models,  wherein  each  crystal  of  a 
polycrystal  is  resolved  geometrically.  Crystal  plasticity 
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models  have  been  implemented  elsewhere  to  study  shock 
loading  of  energetic  materials  cyclotetramethylene  tetranitr- 
amine  (HMX)10  and  pentaerythritol  tetranitrate  (PETN).11'12 
One  aim  of  the  present  work  is  development  and  implemen¬ 
tation  of  a  crystal  plasticity  model  for  RDX,  a  particular 
model  which,  to  the  authors’  knowledge,  has  not  been 
published  elsewhere. 

The  single  crystal  elastic-plastic  model  developed  here 
extends  a  previous  model1  ’  for  cubic  crystals  loaded  to  pos¬ 
sibly  high  pressures.  Here,  the  model  is  applied  to  crystals 
with  orthorhombic  symmetry  characteristic  of  RDX.  Aniso¬ 
tropic  elastic  constants  and  pressure-dependent  compressibil¬ 
ity  are  considered  from  experimental  literature. 14-1 6  Six  slip 
systems  (from  four  different  families  of  systems)  are  imple¬ 
mented  following  analysis  of  indentation  loading  profiles 
and  surface  impressions.6'6  The  model  is  applied  to  study  in¬ 
dentation,  with  a  spherical  indenter  of  (001),  (021),  and 
(210)  faces  of  single  crystals  of  RDX. 

This  paper  is  organized  as  follows:  Constitutive  theory 
and  material  properties  are  described  in  Sec.  II.  Indentation 
simulations  are  reported  in  Sec.  III.  Conclusions  follow  in 
Sec.  IV.  Notation  of  continuum  mechanics  is  used,  e.g.,  bold¬ 
face  type  for  vectors  and  tensors  all  referred  to  fixed  Cartesian 
coordinates.  Summation  applies  over  repeated  indices. 

II.  THEORY 

A.  Single  crystal  model 

Let  x  =  /(X,  t)  denote  the  motion  of  material  points  of 
the  body.  The  deformation  gradient  is 

\7X  =  F  =  FEFP ,  (1) 

where  V(-)  denotes  the  material  gradient  (i.e.,  F„A  =  V  i ya 
=  dxa/dXA),  Fe  denotes  thermoelastic  deformation  of  the 
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crystal  lattice,  and  Fp  represents  plastic  deformation  due  to 
slip.  The  spatial  velocity  gradient  is 

FF~ 1  =  FeFe~  1  +  FeFpFe~  1 ,  (2) 

where  the  plastic  velocity  gradient  associated  with  slip  rates, 
yk,  reference  slip  directions,  Sq,  and  reference  slip  plane  nor¬ 
mals,  ntf],  on  slip  systems,  k,  is 

Lp  =  FpFp-x  =  ^2  yk4  ®  mo-  (3) 

k 


TABLE  I.  Structural  and  physical  properties  of  RDX  single  crystals 
(ambient). 


Property 

Value 

Ref. 

Space  group 

Pbca 

5 

Crystal  structure 

Orthorhombic 

Lattice  parameters  [nm] 

a  =  1.3182 

5 

b=  1.1574 

c=  1.0709 

Mass  density  [g/cm3] 

1.806 

16 

Slip  system  geometry  is  pushed  forward  to  the  current  con¬ 
figuration  via  sk  =  FeSq  and  mk  =  F£~t/Wq,  where  (•) 
denotes  the  transpose.  From  orthonormality  of  slip  directions 
and  slip  plane  normals,  plastic  deformation  is  isochoric; 
thus,  /=  detF  =  det///  >  0.  By  the  polar  decomposition  theo¬ 
rem,  let  Fe  =  Re\Je.  A  logarithmic  thermoelastic  strain  mea¬ 
sure  is  defined  as  £  =  In UE  and  is  split  into  deviatoric  and 
volumetric  parts  as 

£  =  s'  +  (l/3)el,  e  =  tr£  =  ln/,  (4) 

where  1  is  the  second-order  unit  tensor,  and  the  trace  of  a 
second-order  tensor  is  tr(  ).  Let  a  denote  the  usual  Cauchy 
stress  tensor;  stress  in  an  unrotated  coordinate  system  is 
S=RE~laRE. 

Only  the  isothermal  case  is  considered,  an  appropriate 
assumption  for  applications  of  the  model  discussed  in  Sec. 
III.  Let  the  unrotated  stress  be  split  into  deviatoric  and 
hydrostatic  parts, 

S  =  S'  +  S1,  S  =  (l/3)trS  =  (l/3)trff  =  —p,  (5) 

with  p  the  Cauchy  pressure.  The  following  operator  extracts 
the  deviatoric  part  of  a  second-order  tensor, 

F=I  —  (1/3)1  ®  1,  VA  BCD  =  &AC&BD  —  (1/3  )8ab(>CD- 

(6) 

Constitutive  equations  for  deviatoric  stress  and  pressure  are13 

S'  =  l':C:£,+(l/3)(I':C:l)e,  (7) 

p  =  (B0/B') [exp(— B'e)  -  1]  -  (l/3)e'  :  C  :  1.  (8) 

Here,  B0  and  B'  are  the  reference  bulk  modulus  and  pressure 
derivative  of  the  bulk  modulus  and  C  is  the  tensor  of 
second-order  elastic  constants  referred  to  the  (unrotated) 
crystal  frame.  The  colon  denotes  contraction  over  two  pairs 
of  indices,  e.g.,  (C:s')ab  =  Cabcd^ci)  ■  Pressure  dependence 
of  shear  elastic  coefficients,  implemented  elsewhere  for 
cubic  crystals,13  is  omitted  in  Eq.  (8)  because  of  limited  ex¬ 
perimental  data  for  the  material  of  present  interest. 

The  flow  rule  for  slip  is13 

t*  =  *  :  (**  ®  mk)  =  t*[(|t*|  +  £)/y0]msgn(/).  (9) 

The  resolved  shear  stress  on  slip  system  k  is  zk .  Material  pa¬ 
rameters  are  initial  and  constant  slip  strength,  Tq,  for  each 
slip  system,  reference  strain  rate,  y0,  and  rate  sensitivity,  m. 
Constant  ^  <C  y0  provides  a  finite  strength  at  zero  strain  rate. 


B.  RDX 

Physical  properties  of  RDX  single  crystals  are  listed  in 
Table  I.  The  description  applies  to  the  a  phase,  the  stable 
polymorph  for  pressures  under  «3.8  GPa  and  temperatures 
under  «  480K. 

Elastic  properties  are  listed  in  Table  II.  Isentropic 
second-order  elastic  constants15'16  and  bulk  modulus  have 
been  converted  to  isothermal  values  at  295  K  via  the  usual 
thermodynamic  formulae,17  incorporating  anisotropic  ther¬ 
mal  expansion13  and  specific  heat.16  Voigt’s  notation  is 
used:  C abcd  <->  Cap,  where  Greek  indices  1,  2,  ...,  6.  Voigt 
(Gy)  and  Reuss  (G,{)  bounds17  on  the  effective  shear  modu¬ 
lus  are  also  listed.  Differences  between  Voigt  and  Reuss 
bounds  for  bulk  modulus  B0,  on  the  order  1-3%,  are  consid¬ 
ered  insignificant. 

As  is  evident  from  Table  II,  reported  values  of  second- 
order  elastic  constants  can  vary  substantially.  Values 
obtained  using  resonant  ultrasonic  (RUS)  methods15  listed  in 
Table  II  are  in  reasonably  close  agreement  with  those 
reported  by  other  researchers  using  the  same  method.-0  Val¬ 
ues  obtained  using  Brillouin  scattering16  listed  in  Table  II 
are  notably  different,  with  particularly  larger  bulk  stiffness 
and  shear  stiffness  in  certain  directions.  Values  obtained 
using  a  third  technique,  impulsive  stimulated  thermal  scatter¬ 
ing,21  are  similar  to  those  obtained  using  RUS.  Possible  rea¬ 
sons  for  differences  among  measurements  of  elastic 
constants  of  organic  molecular  crystals  are  discussed  else¬ 
where.22  Values  predicted  using  empirical  atomic  mod¬ 
els2’24  also  exhibit  differences  from  those  obtained  in 
experiments,  though  these  predicted  values  tend  to  align 


TABLE  II.  Isothermal  second-order  elastic  constants  of  RDX  (converted 
from  ambient  isentropic  values). 


Property 

Value15 

Value16 

Cu  [GPa] 

24.56 

36.48 

C22  [GPa] 

18.85 

24.49 

C33  [GPa] 

17.33 

20.78 

C12  [GPa] 

7.61 

0.90 

C13  [GPa] 

5.30 

1.26 

C23  [GPa] 

5.24 

8.16 

C44  [GPa] 

5.15 

11.99 

C55  [GPa] 

4.06 

2.72 

C66  [GPa] 

6.90 

7.68 

B0  [GPa] 

10.5 

11.2 

Gv  [GPa] 

6.06 

9.26 

G„  [GPa] 

5.72 

6.40 
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more  closely  with  those  from  RUS  than  from  Brillouin 
scattering. 

In  this  work,  the  two  sets  of  elastic  constants15,16  are 
used,  because  these  are  the  softest  and  stiffest  reported  ex¬ 
perimental  measurements,  respectively;  results  obtained 
using  these  two  sets  might  be  expected  to  bound  the  actual 
response.  Results  obtained  using  these  particular  RUS  con¬ 
stants15  are  very  close  to  those  that  would  be  obtained  using 
similar  values.20  21 

Pressure  and  temperature  dependencies  of  second-order 
elastic  coefficients  have  been  calculated  using  molecular  dy¬ 
namics;2'  however,  these  predicted  values  have  not  been 
validated  using  experiments,  and  some  discrepancies  exist 
among  calculated  second-order  elastic  constants  at  room 
temperature  and  experimental  values.16  For  this  reason,  in 
the  present  study,  the  nonlinear  elastic  model  only  incorpo¬ 
rates  pressure  dependence  of  the  compressibility  obtained 
experimentally 14  and  not  that  of  all  elastic  coefficients,  spe¬ 
cifically,  B'  =  6.95. 

Crystal  plasticity  theories  incorporating  higher-order 
elastic  constants  have  also  been  developed 1_,_;’  and  might 
offer  an  improved  description  of  effects  of  volumetric  and 
shear  deformations  on  tangent  elastic  moduli;  such  an 
approach  is  not  pursued  here,  because  higher-order  elastic 
constants  (e.g.,  20  independent  third-order  constants  for  an 
orthorhombic  crystal17)  are  unknown  for  RDX.  In  a  crystal 
plasticity  model  of  PETN,12  a  tetragonal  crystal,  the  Cauchy 
relations  were  used  to  estimate  unknown  third-order  elastic 
constants. 

Potential  slip  systems  in  RDX — as  identified  from  hard¬ 
ness  versus  orientation  profiles,5  indentation  force  versus 
depth  data,6,7  and  residual  surface  impressions  from  indenta¬ 
tion6 — are  listed  in  Table  III.  Slip  system  geometry  (Fig.  1) 
is  referred  to  a  Cartesian  system  with  axes  (X , ,  X2,  X3)  paral¬ 
lel  to  lattice  vectors  (a,  b,  c ).  Listed  initial  slip  system 
strengths  are  upper  bounds  estimated  from  analysis  of  load 
excursion  data  using  the  analytical  Hertzian  solution  for  fric¬ 
tionless  spherical  indentation  into  a  semi-infinite,  linear  elas¬ 
tic,  isotropic  material.6  Other  systems  may  become  active 
(and  those  listed  may  become  inactive)  for  loading  regimes 
involving  very  different  pressures,  temperatures,  and/or 
strain  rates;  e.g.,  molecular  dynamics  simulations8  suggest 
that  partial  dislocation  loops  may  glide  on  (001)  [010]  during 
shock  loading  at  pressure  in  excess  of  «  1  GPa. 

Strengths  for  various  families  of  systems  are  varied 
parametrically  between  physically  reasonable  bounds  on  the 
order  of  the  theoretical  strength,6,7  where  G  =  GR  =  6.4 


TABLE  III.  Slip  systems  in  RDX  single  crystals  (indentation). 


System  k 

Miller  indices 

m0 

*0 

Max.  strength 
[GPa] 

Ref. 

i 

(021)[100] 

(0,  0.880,  0.475) 

[1,0,  0] 

0.885 

5,6 

2 

(021) [100] 

(0,  -0.880,  0.475) 

[1,0,  0] 

3 

(01 1)[100] 

(0,  0.679,  0.734) 

[1,0,  0] 

0.645 

5,6 

4 

(Oil)  [100] 

(0,  -0.679,  0.734) 

[1,0, 0] 

5 

(010)[100] 

(0,1,0) 

[1,0, 0] 

0.885 

6 

6 

(010)[001] 

(0,1,0) 

[0,  0,  1] 

0.885 

5,6 

{021  }[100]  slip 
2  planes  X  1  direction 


(010)[100]  slip 
1  plane  x  1  direction 


(010)[001]  slip 
1  plane  x  1  direction 


FIG.  1.  Slip  systems  in  RDX  (unit  cell  parameters  not  to  scale). 


GPa.16  Because  strengths  of  families  of  slip  systems  are  not 
known  precisely  a  priori  from  experiments,  shear  strengths, 
Tq,  are  varied  parametrically  over  the  range  listed  in  Table 
IV.  This  range  is  physically  descriptive  of  homogeneous 
nucleation  of  glissile  dislocation  lines  and  loops.6,7,17  Other 
parameters,  which  provide  a  nearly  rate-independent 
response,  are  also  listed  in  Table  IV. 

III.  INDENTATION 
A.  Boundary  value  problem 

The  constitutive  model  of  Sec.  II  is  implemented  in  the 
ALE3D  multi-physics  code.  Simulations  of  indentation  are 
performed  using  an  implicit  solver  for  static  equilibrium. 

The  problem  geometry  mimics  previous  experimental 
studies.6,7  A  spherical  diamond  indenter  of  radius  R  =  1 .482 
pm  is  used  to  indent  a  flat  surface  of  a  single  crystal  of  RDX 
of  variable  lattice  orientation.  Diamond  is  represented  as  an 
isotropic  nonlinear  elastic  material  with  B0  —  443  GPa, 
G  =  538  GPa,  and  B'  =  4.0. 

The  substrate  is  represented  by  a  right  circular  cylin¬ 
der.26  The  cylinder  is  assigned  a  height  and  radius  of  2 R\  fur¬ 
ther  increases  in  dimensions  of  the  cylinder  did  not  affect 
results  of  interest.  The  indenter  is  modeled  as  a  half-sphere. 
Each  body  is  discretized  using  eight-node  hexahedral  ele¬ 
ments  with  selective  reduced  integration.  The  mesh  of  the 
substrate  is  highly  refined  in  the  vicinity  of  contact  beneath 
the  indenter,  where  stress  fields  are  inhomogeneous,  and  the 


TABLE  IV.  Crystal  plasticity  model  parameters  for  RDX. 


Property 

Value 

Tk 

T0 

G/40  <  Tq  <G/10 

m 

5  x  10“3 

ia 

10“2/s 

f 

10_7/s 
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TABLE  V.  Elastic  stiffness  and  Young’s  modulus  of  RDX  oriented  for  in¬ 
dentation  on  planes  (001),  (021),  and  (210). 


Stiffness 

Plane  [RUS15] 

Plane  [Brillouin16] 

(001) 

(021) 

(210) 

(001) 

(021) 

(210) 

Cn  [GPa] 

24.56 

24.56 

20.15 

36.48 

36.48 

22.17 

C22  [GPa] 

18.85 

16.78 

17.33 

24.49 

24.94 

20.78 

C33  [GPa] 

17.33 

17.62 

23.06 

20.78 

26.98 

28.29 

£n  [GPa] 

20.85 

20.85 

16.72 

36.40 

36.40 

19.12 

fi22  [GPa] 

15.69 

14.05 

15.40 

21.28 

22.84 

18.04 

fi33  [GPa] 

15.40 

14.19 

19.31 

18.04 

25.55 

23.61 

mesh  coarsens  progressively  with  distance  from  the  initial 
contact  point.  Simulations  with  smaller  elements  demon¬ 
strated  that  further  increases  in  mesh  refinement  did  not 
affect  results  of  interest. 

During  the  loading  phase,  the  upper  face  of  the  half- 
spherical  indenter  is  assigned  a  constant  (downward)  veloc¬ 
ity  of  D  =  10  nm/s  in  the  laboratory  X3-direction,  leading  to 
strain  rates  similar  to  those  of  experiments.6  The  lower  face 
of  the  cylinder  is  rigidly  fixed,  while  the  lateral  sides  (cir¬ 
cumference)  are  traction-free.  Indentation  depth  is  denoted 
by  D;  actual  depth  d  of  the  tip  of  the  sphere  in  contact  with 
the  surface  is  monitored  as  an  outcome  of  the  solution.  Only 
for  a  rigid  indenter  would  d  =  D.  Indentation  force  P  is  the 
sum  of  nodal  forces  along  the  upper  face  of  the  half  sphere 
acting  in  the  direction  of  D ,  i.e.,  the  sum  of  forces  work  con¬ 
jugate  to  prescribed  nodal  velocities. 

Contact  between  the  indenter  and  substrate  is  assumed 
frictionless,  following  previous  studies  that  rely  on  analytical 
solutions.6  7  Experimental  measurements  of  dynamic  friction 
for  RDX  single  crystals  sliding  on  a  glass  substrate-7  suggest 


a  friction  coefficient  on  the  order  of  unity  for  loads  under  1  g 
(«10  mN),  wherein  contact  is  characterized  as  elastic,  and  a 
friction  coefficient  of  0.35  for  loads  in  excess  of  10  g, 
wherein  contact  is  characterized  as  plastic.  The  present  in¬ 
dentation  simulations  consider  a  different  geometry,  smaller 
system  sizes  (loads  under  1  g),  and  slower  sliding  velocities 
(on  the  order  of  10  nm/s),  so  the  reported  experimental  val¬ 
ues  for  friction  coefficients-7  may  not  strictly  apply  here.  In 
additional  simulations,28  it  was  found  that  differences  in  in¬ 
dentation  force  among  cases  invoking  frictionless  and  stick¬ 
ing  contact  were  insignificant.  In  some  simulations, 
unloading  is  also  performed,  whereby,  after  a  peak  depth  is 
attained,  the  upper  face  of  the  indenter  is  assigned  an  upward 
velocity  of  D=  10  nm/s  until  contact  is  released.  If  plastic 
deformation  has  occurred,  then  some  residual  deformation 
remains  in  the  substrate  upon  unloading. 

The  present  simulations  enable  direct  quantification  of 
surface  and  subsurface  slip  on  each  system.  In  contrast,  hard¬ 
ness  or  indentation  experiments5-7  require  substantial  inter¬ 
pretation  of  data  to  deduce  slip  activity  and  do  not  provide  a 
quantitative  measure  of  relative  contributions  of  each  slip 
system  to  the  overall  strain  field.  In  the  aforementioned 
experiments,  visual  observations  of  slip  traces  are  restricted 
to  residual  surface  profiles,  whereas  simulations  enable  visu¬ 
alization  of  subsurface  slip  activity. 

Simulations  of  indentation  onto  (001),  (021),  and  (210) 
planes  are  reported  in  Sec.  Ill  B.  Let  C<ijk>  be  the  fourth- 
order  matrix  of  elastic  constants  of  the  crystal  oriented  for 
indentation  into  crystallographic  plane  (ijk).  Let  R(llk>  be  the 
corresponding  rotation  matrix.  Then, 


">((/*)  _  n('jk)  ni'jk)  n(ijk)  p(ijk)  ( 

-Alt  Abc  Kr-n  fir  11 


J ABCD 


XAE  KBF  Kcg  KDH  ^EFGH , 


(10) 


FIG.  2.  Indentation  force  vs  applied  dis¬ 
placement;  model  predictions  obtained 
using  RUS  (Ref.  15)  and  Brillouin  (Ref.  16) 
elastic  constants  compared  to  experiment 
(Ref.  6):  (a)  indentation  into  (001);  (b)  in¬ 
dentation  into  (021);  (c)  indentation  into 
(210). 


(c) 
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where  C efgh  corresponds  to  elastic  constants  of  Table  II.  Val¬ 
ues  of  elastic  stiffness  are  listed  in  Table  V.  Perhaps  most  perti¬ 
nent  to  indentation  in  the  ^-direction,  C^10'  >  C^21'1  >  C^01* 
for  either  set  of  constants.  Note  that  from  Brillouin 

scattering16  always  exceeds  the  corresponding  value  from 
RIJ.S. 1  Components  of  compliance  Sxp  are  those  of  the  inverse 
of  Cxp.  Young’s  moduli  are  £11=1/511,  £22=1/522,  and 
£33=1/533  From  Table  V,  from  Brillouin  scattering16 
always  exceeds  the  corresponding  value  from  RUS15  for  a=  1, 
2,  and  3. 

B.  Model  predictions 

Predictions  obtained  using  each  set  of  elastic  con¬ 
stants15,16  are  compared  directly  for  substrate  orientations 
(001),  (021),  and  (210)  in  Fig.  2.  In  results  marked  as 
“elastic”,  slip  is  suppressed  to  permit  assessment  of  the  elas¬ 
ticity  model  in  isolation  and  to  permit  deduction  of  the  yield 
point  upon  comparison  with  “elastic-plastic”  results,  wherein 
the  slip  model  of  Sec.  II  is  enabled.  For  the  latter  results 
shown  in  Fig.  2,  for  each  set  of  elastic  constants,  strength  is 
set  to  a  constant  value  of  Tq  =  G/ 20  for  all  six  slip  systems. 

Significantly  closer  agreement  between  model  and 
experiment6  is  obtained  from  the  elastic  constants  from 
RUS15  for  indentation  onto  (001)  and  (021)  planes.  The  first 
experimental  data  point  in  each  figure  corresponds  to  the 
maximum  depth  at  which  the  indentation  process  remains 
elastically  reversible,  i.e.,  initiation  of  the  first  excursion 
from  a  smooth  force-displacement  profile.6,7  Comparable  ac¬ 
curacy  is  obtained  from  either  set  of  elastic  constants  for  in¬ 
dentation  on  (210)  planes.  The  present  simulations  strongly 
suggest  that  elastic  constants  obtained  from  RUS15,20  pro¬ 
vide  a  more  realistic  representation  of  elastic  stiffness  during 
nano-indentation  than  elastic  constants  obtained  from  Bril¬ 
louin  scattering,16  with  the  latter  appearing  too  stiff.  Results 
shown  in  Fig.  2  invoke  the  geometrically  nonlinear  elastic 
model  with  pressure-dependent  bulk  modulus;  comparison 
with  results  of  additional  calculations28  with  B'  =  0  for  RDX 
demonstrated  that  the  effect  of  nonlinear  compressibility 
becomes  noticeable  for  D>  50  nm. 

In  simulations,  the  yield  point  (initiation  of  slip)  can  be 
deduced  as  the  indentation  depth  beyond  which  elastic  and 
elastic-plastic  model  predictions  begin  to  differ.  From  Fig.  2, 
the  predicted  yield  point  matches  the  experimental  excursion 
point  reasonably  well  for  each  orientation  when  RUS  elastic 
constants15  are  used  in  the  model.  In  contrast,  the  indentation 
depth  at  which  yielding  is  predicted  is  premature  for  (001) 
and  (021)  orientations  when  Brillouin  constants16  are  used. 
This  difference  is  presumably  a  result  of  attainment  of  larger 
resolved  shear  stresses  at  a  given  depth  of  indentation  when 
stiffer  constants  are  used. 

Predicted  forces  exceed  experimental  values  at  larger  in¬ 
dentation  depths  in  each  orientation  and  for  both  sets  of  elas¬ 
tic  constants.  Results  corresponding  to  more  compliant 
elastic  constants16  provide  closer  agreement  to  experimental 
values  than  results  corresponding  to  stiffer  elastic  con¬ 
stants.16  As  noted  elswhere,26  -8  uncertainty  in  the  true  tip  ra¬ 
dius,  R,  of  the  indenter  could  lead  to  discrepancies  between 


simulations  and  experiments.  Surface  fractures  and/or  sub¬ 
surface  fractures  could  contribute  to  a  loss  of  stiffness  that 
would  be  reflected  only  in  the  experimental  data.  RDX  is 
prone  to  cleavage  fracture  on  planes  (001),  (010),  (001), 
(241),  and  (241). 2,4,6  Experimental  data  demonstrate  nearly 
horizontal  steps  in  force  versus  displacement  corresponding 
to  discrete  shear  discontinuities  and/or  fracture  events  that 
are  not  readily  resolved  by  a  constant  strength  continuum 
crystal  plasticity  model,  such  as  the  one  fonnulated  here.  In 


-2.0  -1.0  0.0  1.0  2.0 


(a)  Y  Axis 


(021) 


-2.0  -1.0  0.0  1.0  2.0 
(b)  Y  Axis 


-2.0  -1.0  0.0  1.0  2.0 


(c)  Y  Axis 

FIG.  3.  Cumulative  total  slip  contours  for  indentation  to  depth  D  —  200  nm 
using  elastic  constants  from  RUS  (Ref.  15)  and  slip  system  strength  G/ 20;  a 
slice  along  the  centerline  of  the  cylinder  normal  to  the  laboratory  X\  axis  is 
shown:  (a)  indentation  into  (001)  [X-axis  normal  to  (100)];  (b)  indentation 
into  (021)  [X-axis  normal  to  (100)];  (c)  indentation  into  (2 10) [X-axis  normal 
to  (230)]. 
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TABLE  VI.  Maximum  local  slip  yk  at  indentation  depth  of  200  nm  for  in¬ 
dentation  on  (001),  (021),  and  (210)  planes;  RUS  elastic  constants.15 


System  k 

Plane  (001) 

Plane  (021) 

Plane  (210) 

i 

0.01 

0.19 

0.40 

2 

0.01 

0.00 

0.40 

3 

0.08 

0.16 

0.29 

4 

0.08 

0.00 

0.29 

5 

0.00 

0.10 

0.51 

6 

0.31 

0.41 

0.16 

all  (y) 

0.31 

0.46 

0.52 

experiments,6  surface  fractures  are  detected  only  at  loads  sig¬ 
nificantly  exceeding  yield  excursion. 

Details  of  slip  system  interactions  and  pressure  depend¬ 
ence  of  shear  strength  are  omitted  in  the  present  model.  In¬ 
dentation  experiments  have  suggested  the  importance  of 
cross  slip,6  and  atomic  modeling  has  noted  that  different  slip 
mechanisms  may  become  important  at  high  pressures.8 
Incorporation  of  these  effects  into  a  more  complex  slip 
model  might  provide  closer  agreement  with  experiment,  e.g., 
if  glide  resistance  were  to  decrease  with  pressure,  since  local 
pressures  under  the  indenter  can  achieve  several  GPa.28 
Atomic  modeling8'24  may  provide  insight  into  dependence  of 
slip  resistance  on  pressure  (and  temperature,  etc.)  not  avail¬ 
able  from  experimental  methods. 

Model  predictions  of  cumulative  slip  for  indentation  to  a 
depth  of  D  =  200  nm  onto  planes  (001),  (021),  and  (210)  are 
shown  in  Fig.  3.  In  each  case,  a  slice  along  the  centerline  of 
the  cylinder  normal  to  the  laboratory  Xx  axis  (i.e.,  X-axis)  is 
shown.  For  indentation  on  (001)  and  (021)  planes,  the  X , 
axis  is  normal  to  a  (100)  plane;  for  indentation  on  (210),  the 
Xi  axis  is  approximately  normal  to  a  (230)  plane.  Cumula¬ 
tive  total  slip,  y,  is  defined  as 

y  =  =  (n) 

k  k  •’ 

with  yk  the  monotonically  increasing  cumulative  slip  on  sys¬ 
tem  k.  Slip  activity  is  greater  for  indentation  into  (021)  and 
(210)  planes  than  for  indentation  on  (001).  For  this  particular 
viewing  plane,  slip  contours  for  indentation  on  (021)  are 
noticeably  asymmetric.  The  wireframe  mesh  of  the  indenter 
is  drawn  in  each  figure;  in  order  to  enable  clear  visualization 
of  heterogeneous  slip  distributions  in  the  RDX,  the  mesh  of 


TABLE  VII.  Maximum  local  slip  yk  at  indentation  depth  of  200  nm  for  in¬ 
dentation  on  (001),  (021),  and  (210)  planes;  BriUouin  elastic  constants.16 


System  k 

Plane  (001) 

Plane  (021) 

Plane  (210) 

i 

0.01 

0.19 

0.41 

2 

0.01 

0.00 

0.41 

3 

0.06 

0.13 

0.28 

4 

0.06 

0.00 

0.28 

5 

0.01 

0.10 

0.52 

6 

0.35 

0.51 

0.26 

all  (y) 

0.35 

0.56 

0.58 

the  substrate,  which  is  considerably  more  refined  than  that  of 
the  indenter,  is  not  shown. 

Tables  VI  and  VII  list  maximum  local  cumulative  slip 
(i.e.,  the  maximum  value  of  yk  at  any  location  in  the  RDX 
substrate)  at  an  indentation  depth  D  =  200  nm  when  various 
elastic  constants1516  are  implemented.  Total  slip,  y,  listed  in 
the  bottom  row  of  each  table  is  not  necessarily  the  sum  of  all 
yk  listed  in  a  given  column,  because  the  location  in  the  sub¬ 
strate  where  total  slip  is  maximum  does  not  necessarily  cor¬ 
respond  to  the  location  where  each  yk  is  maximum. 
Orientation  (001)  exhibits  slip  primarily  on  system  6.  Signifi¬ 
cant  activity  of  multiple  slip  systems  is  evident  for  indenta¬ 
tion  onto  (021)  and  (210)  planes.  Trends  are  qualitatively 
similar,  regardless  of  choice  of  elastic  constants.  Cumulative 
slip  magnitudes  are  generally  slightly  larger  in  Table  VII 
than  in  Table  VI,  because  larger  stresses  are  attained  at  the 
same  depth  of  indentation  when  stiffer  elastic  constants  are 
prescribed. 

Elastic-plastic  simulations  were  also  performed,  wherein 
the  strength,  Tq,  of  one  family  of  systems  was  set  to  G/ 20, 
with  that  of  all  others  set  to  G/10.  Comparison  of  force  ver¬ 
sus  depth  profiles  among  these  simulations  enabled  further 
assessment  of  the  most  active  slip  systems  for  each  crystal 
orientation.28  Results  were  consistent  with  those  in  Tables 
VI  and  VII:  slip  system  6  is  dominant  for  indentation  on 
(001)  and  (021),  while  systems  1-5  are  all  important  for  in¬ 
dentation  on  (210). 

Model  predictions  for  uniform  slip  strengths,  Tq,  of  G/ 
10,  G/ 20,  and  G/40  are  compared  with  experimental  data6  in 
Fig.  4  for  indentation  on  (210).  Reduction  of  slip  strength 
from  G/ 20  to  G/40  enables  closer  agreement  with  experiment 
at  larger  indentation  depths,  but  also  leads  to  premature  ini¬ 
tiation  of  slip  and  under-prediction  of  force  at  smaller  depths 


FIG.  4.  Indentation  force  vs  applied  dis¬ 
placement  for  elastic-plastic  indentation  into 
plane  (210)  with  uniform  slip  strengths  G/ 
10,  G/ 20,  and  G/40  compared  to  experiment 
(Ref.  6);  model  predictions  obtained  using 
(a)  RUS  (Ref.  15)  and  (b)  Brillouin  (Ref.  16) 
elastic  constants. 


(a) 


(b) 
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FIG.  5.  Indentation  force  vs  applied  displacement  for  a  single  load-unload 
cycle  on  each  of  (001),  (021),  and  (210)  planes;  elastic  constants  from  RUS 
(Ref.  15),  uniform  slip  strength  G/ 20. 

near  experimental  excursion  points,  i.e.,  near  initial  yield. 
Similar  trends  were  predicted28  for  indentation  onto  (001) 
and  (021)  planes  when  strength  is  increased  or  decreased 
uniformly  among  slip  systems. 

Loading-unloading  simulations  of  indentation  on 
(001),  (021),  and  (210)  planes  were  performed  to  enable 
comparison  with  experimental  observations.  Force  versus 
displacement  predictions  are  shown  in  Fig.  5  for  a  single 
load-unload  cycle,  corresponding  to  loading  to  D  «  200  nm 
and  then  unloading  to  D  =  0.  During  much  of  the  loading 
phase,  orientation  (001)  is  slightly  less  stiff  than  orienta¬ 
tions  (021)  and  (210).  However,  at  Dm  200  nm,  P  is 
largest  for  orientation  (001).  Recall  that  orientation  (001) 
demonstrates  a  lower  elastic  stiffness,  but  also  less  slip 


activity.  Hence,  at  larger  indentation  depths,  increased  slip 
activity  for  orientations  (021)  and  (210)  lowers  their 
effective  tangent  stiffness  and,  hence,  P  below  that  of 
orientation  (001).  Hysteresis  is  also  substantially  greater 
for  indentation  on  (021)  and  (210)  than  for  (001),  again 
demonstrating  less  slip  activity  in  the  latter.  Orientation 
(210)  demonstrates  the  most  hysteresis  (and  greatest  slip 
activity)  and  the  largest  elastic  stiffness  for  much  of  the 
unloading  phase. 

Figure  6  shows  total  cumulative  slip,  y,  on  unloaded 
surfaces  of  (001),  (021),  and  (210)  planes.  In  each  case,  the 
maximum  indentation  depth  prior  to  unloading  is  D  «  200 
nm.  Non-circular  contours  are  consistent  with  activity  of 
fewer  than  five  geometrically  independent  systems  necessary 
to  accommodate  an  arbitrary  plastic  strain  field.6 

Surface  slip  contours  from  individual  systems  were  also 
examined.28  Predicted  surface  slip  activity  results  primarily 
from  system  6  for  indentation  on  (001).  Slight  contributions 
to  the  circular  total  slip  trace  in  Fig.  6(a)  are  due  to  systems 
3  and  4,  i.e.,  {01 1 } [  100] .  Surface  slip  activity  is  predomi¬ 
nantly  from  system  6,  with  minor  contributions  from  systems 
1  and  3,  for  indentation  on  (021).  Surface  slip  activity  is  pre¬ 
dominantly  from  system  5,  with  contributions  from  systems 
1,  2,  and  6,  for  indentation  on  (210).  Faint  contributions 
from  systems  3  and  4  are  also  predicted  for  indentation  on 
(210). 

Active  slip  planes  during  the  loading  history  at  the  spec¬ 
imen  surface,  as  predicted  here,  are  compared  with  those 
deduced  from  experimental  surface  impressions6  in  Table 
VIII.  It  is  noted  that,  in  simulations  of  indentation  on  (210), 


FIG.  6.  Residual  total  slip  contours  at  the 
surface  for  indentation  to  depth  D  ~  200  nm 
followed  by  unloading;  results  obtained 
using  elastic  constants  from  RUS  (Ref.  15) 
and  uniform  slip  system  strength  G/ 20:  (a) 
indentation  into  (001);  (b)  indentation  into 
(021);  (c)  indentation  into  (210). 


(c) 
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TABLE  VIII.  Significantly  active  slip  planes  at  the  specimen  surface  during 
indentation  of  (001),  (021),  and  (210)  planes. 


Indentation  surface 

Present  results 

Experiment6 

(001) 

(010),  {011} 

(010),  {Oil},  {021} 

(021) 

(010),  (011),  (021) 

Not  reported 

(210) 

(010),  {Oil},  {021} 

(010),  {Oil},  (021) 

systems  3  and  4,  i.e.,  {011}  [100] ,  demonstrate  substantial 
activity  within  the  bulk  material  (i.e.,  beneath  the  surface), 
but  little  cumulative  slip  at  the  surface.  Though  not  presented 
graphically  in  this  paper,  predicted  residual  bulk  and  surface 
slip  contours  were  also  obtained  using  the  relatively  stiff 
elastic  constants  from  Brillouin  scattering.16  Predicted  slip 
contours  appeared  qualitatively  similar,  regardless  of  choice 
of  elastic  constants  from  RUS15  or  Brillouin.16 


IV.  CONCLUSIONS 

A  nonlinear  anisotropic  elastic-plastic  model  has  been 
developed  for  RDX.  The  model  accounts  for  orthorhombic 
elastic  constants,  pressure-dependent  compressibility,  and 
dislocation  glide  on  up  to  six  distinct  slip  systems. 

Numerical  simulations  of  spherical  indentation  on 
(001),  (021),  and  (210)  planes  of  single  crystals  show  sig¬ 
nificant  influences  of  elastic  anisotropy  and  nonlinearity  on 
force-displacement  data.  Model  predictions  for  initial  elas¬ 
tic  response  using  constants  measured  with  resonant  ultra¬ 
sound  spectroscopy  agree  with  experimental  force- 
displacement  data  for  indentation  on  (001),  (021),  and 
(210)  planes.  Predictions  using  constants  measured  with 
Brillouin  scattering  are  in  reasonable  agreement  with 
experiments  for  indentation  on  (210),  but  are  stiffer  than 
experiments  for  indentation  on  (001)  and  (021).  Orientation 
(001)  is  elastically  most  compliant,  in  agreement  with 
experiments. 

Critical  shear  strengths  associated  with  slip  initiation 
have  been  estimated  as  G/ 20,  where  G  is  a  representative 
elastic  shear  modulus.  Initial  yield  points  predicted  by  the 
model  are  in  close  agreement  with  experimental  load  excur¬ 
sion  data  when  elastic  constants  from  resonant  ultrasound 
spectroscopy  are  used.  Predictions  of  force  for  larger  inden¬ 
tation  depths,  wherein  predicted  plastic  slip  is  substantial, 
tend  to  exceed  experimental  values,  regardless  of  which  set 
of  elastic  constants  is  used.  The  constant  strength  (i.e.,  per¬ 
fectly  plastic)  slip  model  implemented  here  is  unable  to  rep¬ 
licate  nearly  horizontal  steps  in  indentation  force  observed  in 
experiments.  Such  steps  may  correspond  to  discrete  slip 
events  of  width  too  fine  to  be  captured  by  a  conventional 
continuum  slip  model.  Dependencies  of  shear  strength  on 
slip  history  and  pressure  have  been  omitted;  incorporation  of 
such  physics,  for  example,  as  suggested  by  atomic  simula¬ 
tions,  might  provide  improved  agreement.  Fractures 
observed  in  experiments  (at  the  surface)  or  not  observed 
(subsurface)  would  also  explain  the  higher  stiffness  in  simu¬ 
lations  relative  to  experiments. 

Simulations  suggest  that  slip  planes  (010)  and  {011} 
contain  active  systems  for  indentation  on  (001),  with  slip  on 


system  (010)[001]  dominating  the  inelastic  response;  experi¬ 
mental  surface  observations  confirm  that  these,  as  well  as 
{021}  slip  planes,  may  also  be  active.  Simulations  suggest 
that  slip  planes  (010)  and  {021},  and  to  a  lesser  extent 
{Oil},  are  active  at  the  specimen  surface  for  indentation  on 
(210);  these  same  planes  have  also  been  confirmed  as  active 
in  experiments.  Simulations  suggest  that  planes  (010),  (01 1), 
and  (021)  contain  active  systems  for  indentation  on  (021); 
particular  slip  planes  active  for  this  orientation  have  not  been 
reported  from  experiments. 

The  present  results  suggest  that  system  (010)[001  J  pro¬ 
vides  the  largest  contribution  to  the  inelastic  material 
response  (i.e.,  post-yield  force  versus  displacement  curve) 
for  indentation  on  (001)  and  (021)  planes,  while  five  systems 
{021} [100],  {01 1}[100],  and  (010)[100]  all  contribute  to 
inelastic  response  for  indentation  on  (210)  planes.  Plastic  de¬ 
formation  and  hysteresis  are  more  extensive  for  indentation 
on  (021)  and  (210)  than  for  indentation  on  (001).  Since  much 
plastic  deformation  occurs  in  the  bulk  of  the  material,  and 
since  different  slip  mechanisms  may  be  prominent  at  the  sur¬ 
face  and  in  the  bulk,  the  present  results  offer  new  insight  into 
inelastic  mechanical  behavior  of  RDX  not  available  from  ex¬ 
perimental  observations  of  residual  surface  topography 
alone. 

The  model  developed  here,  when  used  with  elastic  con¬ 
stants  obtained  from  resonant  ultrasonic  methods,  is 
thought  to  provide  an  accurate  representation  of  the  nonlin¬ 
ear  anisotropic  response  of  RDX  single  crystals  up  to  and 
including  the  onset  of  slip.  The  present  model  is  also 
thought  to  provide  a  qualitatively  reasonable  depiction  of 
activity  of  different  slip  systems  when  a  uniform  and  con¬ 
stant  shear  strength  on  the  order  of  G/ 20  is  prescribed. 
Refinements  of  the  model  are  needed  to  address  any  reduc¬ 
tion  in  stiffness  associated  with  discrete  or  highly  localized 
shear  events  or  cleavage  fractures  observed  at  larger  inden¬ 
tation  depths. 
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